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ABSTRACT 

In the decade of the dawn of gravitational wave astronomy, the concept of multimes- 
senger astronomy, combining gravitational wave signals to conventional electromag- 
netic observation, has attracted the attention of the astrophysical community. So far, 
most of the effort has been focused on ground and space based laser interferometer 
sources, with little attention devoted to the ongoing and upcoming pulsar timing arrays 
(PTAs). We argue in this paper that PTA sources, being very massive (> 1O 8 M0), 
cosmologically nearby (z < 1) black hole binaries (MBHBs), are particularly appealing 
multimessenger carriers. According to current models for massive black hole formation 
and evolution, the planned Square Kilometer Array (SKA) will observe thousands of 
such massive systems, being able to individually resolve and locate in the sky several 
of them (maybe up to a hundred). MBHBs form in galaxy mergers, which are usu- 
ally accompanied by strong inflows of gas in the center of the merger remnant. By 
employing a standard model for the evolution of MBHBs in circumbinary discs, with 
the aid of dedicated numerical simulations, we characterize the gas-binary interplay, 
identifying possible electromagnetic signatures of the PTA sources. We concentrate 
our investigation on two particularly promising scenarios in the high energy domain, 
namely, the detection of X-ray periodic variability and of double broad Ka iron lines. 
Up to several hundreds of periodic X-ray sources with a flux > 10 _13 erg s _1 cm -2 will 
be in the reach of upcoming X-ray observatories; in the most optimistic case, few of 
them may be already being observed by the MAXI detector placed on the Interna- 
tional Space Station. Double relativistic Ka lines may be observable in a handful of 
low redshift (z < 0.3) sources by proposed deep X-ray probes, such as Athena. The 
exact figures depend on the details of the adopted MBHB population and on the prop- 
erties of the circumbinary discs, but the existence of a sizable population of sources 
suitable to multimessenger astronomy is a robust prediction of our investigation. 

Key words: Black hole physics - Accretion, accretion discs - Gravitational waves - 
pulsars: general - X-rays: general 



1 INTRODUCTION 

Within this decade the detection of gravitational waves 
(GWs) may be a reality, opening a completely new win- 
dow on the Universe. While signals coming from compact 
stars and binaries fall in the observational domain of op- 
erating and planned ground based interferometers (such as 
LIGO, VIRGO, and the proposed Einstein Telescope), mas- 
sive black hole (MBH) binaries (MBHBs) are expected to 
be among the primary actors on the upcoming low fre- 
quency stage, where the 10 -4 — 10 _1 Hz window is going 
to be probed by gravitational wave interferometry in space 



(see, e.g., the proposed Laser Interferometer Space Antenna 
(LISA)). Moving to even lower frequencies, long term moni- 
toring of an array of millisecond pulsars (forming a so called 
pulsar timing array, PTA) may unveil the characteristic fin- 
gerpr int left by GWs in the time of arrival of the radio pulses 
(e.g.. lSazhinll978l : lHellings fc DownJl983h. T he Parkes Pul- 
sar Timing Array fPPTA. lManchesteiJl2008lb the E uropean 
Pulsar Timing Array fEPTA, Ijanssen et aUl200Sf ) and the 
North American Nan ohertz Obse rvatory for Gravitational 
Waves (NANOGrav, I Jenetl 120091 ). joining t ogether in th e 
International Pulsar Timing Array (IPTA, iHobbsllioiol '), 
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are already collecting data and improving their sensitivity 
in the frequency range of ~ 10 -9 — 10 -6 Hz, and in th e next 
decad e the planned Square Kilometer Array (SKA, lLaziol 
120091 ) will provide a major leap in sensitivity. Even though 
GW observations alone will be an outstanding breakthrough 
in science, the astrophysical payouts will be greatly ampli- 
fied by the coincident identification of electromagnetic coun- 
terparts. The identification of the host galaxy of a MBH bi- 
nary merger will (i) improve our understanding of the nature 
of the galaxy hosting coalescing MBHBs (e.g. galaxy type, 
colours, morphology, etc.); (ii) help to reconstruct the dy- 
namics of the merging galaxies and of their MBHs; and (iii) 
offer the possibility of studying accretion phenomena onto 
systems of known mass and spin (which, also in PT A obser- 
vations, can be measured from the GW signal, e.g. | Vecchiol 
12004 ISesana fc Vecchidl20ld : ICorbin fc Cornishll2010l) ~ 

Electromagnetic counterparts to GW events have at- 
tracted a lot of attention in the last few years in the context 
of LI SA observations of coalescing MBHBs (see ISchnittmanl 
l201ll . and references therein). The final coalescence being a 
violent event, many possible counterparts related to shocks 
flares and transients in the surrounding ambient gas have 
been proposed. However, the bulk of the LISA sources, are 
likely to cluster at 3 < z < 8 f Ses ana et al.ll2007T l. Given the 
modest s ky localization accurac y of the detector (generally 
> ldeg 2 , lLang fc Hughes! I2008T ). the number of candidate 
hosts in the putative GW errorbox is likely to be very high 
(order of millions). Moreover, our incomplete understanding 
of the observational signatures of such events, together with 
their intrinsic weakness (most LISA sources are expect ed to 
be MBHBs with masses < 10 6 M Q ISesana et all 120071 ) will 
make host identification extremely challenging. 

On the other hand, very little attention has been de- 
voted to possible counterparts of PTA sources. An obvi- 
ous advantage with respect to LISA sources is that those 
are expected to be very massive (M > 10 8 Mp)), c osmo- 
logically nearby (z < 1) systems l|Sesana et al.l 12009). Any 
characteristic electromagnetic signature would be there- 
fore within the sensitivity range of current and future ob- 
servation capabilities. PTA sources are inspiralling MB- 
HBs still far from coalescence (emitting in the ~ 10~ 9 — 
10~ 6 Hz frequency range), that can be considered, for 
any practical purpose, sta tionary during typical observa- 
tion timescales of decades IjSesana fc Vecchioll2010P l. More- 
over, fairly high eccentricity should be the no rm, as shown 
by MBHB harden i ng studies in both stellar (ISesana 20101 ; 



ISesana et al.ll201l];lKhan et alJl201ll:|Preto et alj|2011) and 



gaseous JArmitage fc ^^^H S ^^Tll 



2009; 



iRoedig et al.l l201ll ) environments. If the binary is sur- 
rounded by a circ umbinary gaseous disc, eccen t ricity triggers 
perio dic inflows IjArtvmowicz fc Lubowill996l ; IRoedig et al.l 
l201ll ). opening a wide range of possibility for distinctive bi- 
nary fingerprints that can be observationally identified. 



I2005T ) databasqj, we construct detailed MBHB populations 
satisfying all the currently available constraint s in terms of 
local MB H mass function (iMarconi et alj|2004f ). MBH-host 
relations l|Gultekin et al. 2009), and MBHB m erger rates 
as in ferred from close galaxy pair counts (e.g. iBell et al.l 
120061 ). We assume that after a galaxy merger, cold gas is 
funneled in t he center of the remnant , forming a circum - 
nuclear disc (JMihos fc Hernquistl [l996l ; iMaver et ai] 120071 ). 
The new-formed sub-parsec MBHB excavates a hollow re- 
gion in the center of the disc (usually referred to as 'gap'), 
and the dynamics of the system is governed by the mutual 
MBHB-disc torques. We compute the detailed population of 
GW emitting sources at each frequency by adopting a simple 
analytical model for their evolution under the dynamical ef- 
fect o f gaseous drag and GW emission fsee lKocsis fc Sesanal 
l201ll . for details). After identifying the prototypical MBHB 
observable with PTAs, we simulate its dynamics numerically 
using a modified version of the smooth particle hydrodynam- 
ics (S PH) code GADGET-2 (|Springell [JOOJI ; ICuadra et all 
I2009T ) . to track the fate of the material leaking through the 
gap and captured by the MBHs. The outcome of the simula- 
tion are then used to model analytically the accretion onto 
the two MBHs and the emitted radiation. 

Several distinctive signatures of inspiralling MBHBs 
have been presented in the literature. Proposed scenar- 
ios range from the radio and the optical, up to the 
X-ray s, including variabi lity of the optical continuum 
(e.g., iHaiman et al.l |2009|). spectra l shifts in the broad 
line emission (IBegelman et al.l 1 19801 ; lEracleous et al.l l201ll ; 
iTsalmantza et al.ll201ll ). pec uliar flux ratios betw een opti- 
cal/UV broad emissio n lines flM ontuori et alj|201ll), precess - 
ing or x-shaped je ts dLiul [20041; lLobanov fc Roland! 120051 ). 



periodic outbursts IjSillanpaa et al.l 19881). astrom etric mea- 



surement of the source motion i Sudou et al.ll2003r i. Yet, the 
only unambiguous candidate is a double compac t radio core 
at a 7pc projected separation JRodriguez et al.l IJ2006T ). We 
identify here a number of possible characteristic signatures 
and we concentrate our investigation on two particularly 
promising scenarios in the high energy domain, namely, the 
detection of X-ray periodic variability and of double broad 
Kq iron lines. For both scenarios, we quantify the popula- 
tion of potentially observable sources, and the joint detection 
prospects with future X-ray observatories and PTAs. 

The paper is organized as follows. In Section 2 and 3 we 
model and describe in detail the MBHB population relevant 
to PTA observations; whereas in Section 4 we investigate the 
dynamics of the typical source by combining high resolution 
SPH simulations to an analytical model for the accreted 
matter. In Section 5 we model the emitted spectrum, iden- 
tifying possible characteristic signatures, and Section 6 and 
7 are devoted to periodic variability and double broad Ka 
lines, respectively. We discuss strategy for joint PTA and 
X-ray detection and draw our conclusions in Section 8. 



In this paper we study the prospects of conducting 
multimessenger astronomy with PTA observations of MB- 
HBs. Our aim is to quantify the number of sources possi- 
bly detectable both in the GW and in the electromagnetic 
realm, identifying their signatures, to forecas t future identi- 
fication. Starting from the Millennium Run (|Springel et al.l 



1 http://www.mpa-garching.mpg.de/galform/virgo/millennium/ 
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2 THE PULSAR TIMING MASSIVE BLACK 
HOLE BINARY POPULATION 

Our main goal is to quantify and characterize the pop- 
ulation of MBHBs that is accessible to observation both 
in GWs through pulsar timing and in the electromagnetic 
realm via distinctive signatures. Formally, what we need 
to estimate is the population of MBHBs observable to- 
day in the universe as a function of MBH masses, redshift 
and orbital frequency (or, alternatively, semimajor axis): 
d 4 N/(dMidM2dzd\nf r , k ). Here M\ > M 2 are the masses 
of the two binary components, z is the source redshift 
and f r ,k is the rest-frame Keplerian frequency of the sys- 
tem. This can be formally written as the product of (i) 
d 4 JV/ (dMidM^dzdtr) , times (ii) di r /dln/,..fc. Item (i) is noth- 
ing else but the cosmological coalescence rate of MBHBs, 
and it depends on the general clustering history of struc- 
tures in the Universe; item (ii) describes the time evolution 
of the binary frequency, and is determined by the detailed 
physical processes driving the binary dynamics prior to coa- 
lescenceQ Once such a population has been determined, the 
requirement of observability through PTA, selects a subset 
of systems which we will refer to as the PTA-MBHB popu- 
lation. In the following we will treat items (i)-(ii) separately 
(Sections 12.11 and 12. 2[) . and we will describe the character- 
istic timing residuals induced in a PTA, defining the PTA- 
MBHB population (Section [23)) . 



2.1 Cosmological coalescence rate 

To construct the MBHB cosmological coalesc ence rate, we 
follow the procedure described in Section 2 of ISesana et al.l 
(2009), the reader is deferred to that paper for full details. 
In few words, we extract cat alogs of mergin g galax ies from 
the semi-analytical mod el of Bertonc et al. 1J2007I ) applied 
to the Millennium run (jSpringel et al.l 120051 ). We then as- 
sociate a central MBH to each merging galaxy in our cata- 
log ue. We adopt t h e rece nt fit to the M — a relation given 
by iGiiltekin et al.l (J2009) , and we consider accretion to be 
efficient onto both MBHs before the final coalescence of the 
binary. 

Assigning a MBH to each galaxy, we obtain a catalogue 
of all the mergers occurring in the (500/ftMpc) 3 comoving 
volume of the simulation, labeled by MBH masses and red- 
shift. From this, we numerically generate the cosmic merger 



rate, d N / ' (dM\dMidzdtr) , by weighting each event with the 
observable comoving volume shell at every redshift. 



2.2 MBH binaries in circumbinary discs 

We assume that MBHBs ev olve in geom e trical ly thin 
circumbinary accretion discs IjHaiman et al.l 120091 ). Even 
though such an assumption is enforced by the requirement 
of a distinctive electromagnetic signature, we find that most 



2 For a given coalescence rate, dt r /dlnf r j, quantifies the number 
of binaries as a function of orbital frequency that has to be present 
at any time in the sky to sustain that particular rate. This item 
can be estimated taking into account that we are interested in 
subparscc MBHB, for which we can employ simple evolutionary 
models driven by disc-binary mutual torques and GW emission, 
as we will see in Section 2.2. 



of the merging galaxy pairs extracted by the Millennium run 
involve at least one massive, gas rich, spiral galaxy. Signifi- 
cant in flows of cold gas are therefo re expected in the remnant 
nuclei (jMihos fe Hernquistlll996t ). feeding the circumbinary 
disc. 

Before entering the dynamics in detail, we define all the 
quantities used for describing the MBHB-disc system. The 
MBHB has total mass M = Mi + M 2 (Mi > M 2 ), mass 
ratio q = M1/M2, symmetric binary mass ratio q s — 4g/(l-|- 
q 2 ), semimajor axis a, rest-frame orbital Keplerian frequency 
fk,r, and eccentricity e. The thin disc is described in terms 
of the a viscosity parameter, the accretion rate M and the 
accreted mat t er/rad iation conversion efficiency e. Following 
iRoedig et al.l (|2011f ). we define 5 to be the relative size of 
the gap to the binary semimajor axis, and we take a fiducial 
value S — 2. Throughout the paper, expressions are given 
in units of M 8 = M/10 S M Q , a . 3 = a/0.3, e .i = e/0.1 
and rho.3 = m/0.3, where rh = M /M-Edd is the accretion 
rate normalized to the Eddington rate. Subscript '3' refers 
to lengths given in units of 10 3 Rs where R s = 2GM/c 2 is 
the Schwarzschild radius associated to the total mass of the 
binary. 

The picture we adopt for the MBHB evolution 
is the following. At separations relevant to our work 
(a < 0.0 3 pc), the binary has excav ated a cavity (gap, 
see, e.g. lArtvmowicz fc Lubowl [1994) in the circumbi- 
nary acc retion disc. Both viscous torques exerted by 
the disc ([G oldrcic h fc Tremaind [1980) and GW emission 
IjPeters fc Mathews! 1 19631 ) dissipate the binary binding en- 
ergy causing its orbital decay. The secondary MBH is, in 
general, more massive than the local disc and the problem 
is analogous to secondary-dominated Type-II migration in 
planetary dynamics. A self-consistent solut ion to this prob- 
lem was provided by ISver fc Clarke! |l995j) , the interested 
reader is referred to their original paper and to lHaiman et al.l 
(2009) and Kocsis & Scsana (2011) for a detailed discussion 
in the context of MBHB evolution. Such a solution applies 
only to discs with decreasing s urface density as a functio n 
of the radius. Standard a-discs (jShakura fc SunvaevNl973r ). 
with viscosity proportional to the total pressure, meet such 
a requirement only in the gas pressure dominated zone, but 
not in the inner radiation pressure dominated zone. /3-discs, 
with viscosity proportional to the gas pressure only, fulfill 
this condition at all radii. We therefore assume /3-discs as 
our fiducial disc model for describing the binary-disc mi- 
gration. All the relev ant features of /3-discs are reviewed by 
iHaiman et al.l (|2009|), here we merely introduce the quan- 
tities relevant to the practical computation of the PTA- 
MBHB population. The relevant timescales in the system 
are: 
(i) the viscous timescale 

*„ = 6.09 x 10 5 yra-j /5 ( *M )"'* M 8 6/5 5 7/5 a 3 7/5 ; (1) 



eo.i 



-5/8 



(ii) the migration timescale 



; nnn, m6 -1/2 / m 0.3 

t m — 2.09 x 10 yra ,' 

eo.i 



(iii) the GW shrinking timescale 
few = 7.84 x 10 7 yrM 8 g s " 1 a3F , (e)" 1 
where 



(3) 
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F(e) = (l-e 2 ) 



-7/2 



, 73 2 37 4 

1 H e H e 

24 96 



(4) 



These timescales define two characteristic binary separa- 
tions. Equating t m to i gw we get the decoupling separation 



dec n o-i -4/25 / "T-0.3 

a 3 =0.31a .3' 

£0.1 



-1/5 



M -V25 g ll/25 ( jr/25 F(e) 8/26 ) (5) 



which is the separation below which the GW driven migra- 
tion is faster than the gas driven migration, and the binary 
evolution is dictated by GW emission only. This is the rele- 
vant separation that has to be considered when computing 
the frequency distribution of inspiralling MBHBs. Equating 
t v to t gw we get the 'disc freezing' (or detachment) separa- 
tion 



fr n on -4/13 

a 3 =0.22a 03 / 



m .3 
eo.i 



-2/13 



M 8 1/13 g s 4/13 5 7/13 F(e) 5/13 ,(6) 



below which the binary shrinking timescale is faster than the 
viscous inward diffusion of the inner edge of the disc. After 
this point, the disc can no longer follow the binary, and 
it is basically 'frozen' during the quick subsequent inspiral 
leading to the MBHB coalescence. This latter characteristic 
radius discriminates between binaries attached to their discs 
and binaries detached from their discs, and it has to be con- 
sidered in computing the population of 'observable systems'. 
In general a T < a dcc , defining three different phases relevant 
to our study: 

• phasel, a > a dcc . The binary is coupled to the disc 
which regulates its dynamical evolution; 

• phasell, a iT < a < a dcc . The binary is dynamically de- 
coupled from the disc and it is driven by GW emission. How- 
ever, the viscous time is short enough that the disc can follow 
the binary in its inspiral; 

• phaselll, a < a iT . The binary leaves the disc behind and 
quickly coalesces. 

During phasel, th e binary maintains a constant limiting 
eccentricity given by (|Roedig et al.ll2011f ) 



eo 



0.66Vln5- 0.65 + 0.19 



(7) 



while, after decoupling, in phasell and III, the eccentricity 
can be numerically computed by solving the implicit equa- 
tion 



/k,r 
./dec 



1 



e 

"dec 



l + ±2±e 2 

' 304 



1+ i2ie 2 

T 304 dec 



-STO ' 



-3/2 



(8) 



obtained in the quadrupole approximation by coupling 
the GW orbital decay r ate to the eccentricity decay rate 
l|Peters fc Mathews! fl963h . 

To compute t he frequency dist r ibutio n of inspiralling 
MBHBs, following iKocsis fc Sesana! (J201l|), we express the 
frequency evolution of the binary in terms of 'residence time' 



lit r 



dlnfrj 



dt r dlna 



dlna dln/ r . 



_ 2 



(9) 



In phasel, when the binary is driven by the disc, i res = 
t m ; while in phasell and phaselll, GW emission takes over 
and ires = tgw. Putting the pieces together, for any given 
MBHB-disc parameters, we compute a dcc , assuming eo given 
by equation ([7]); dt r j d\nf r ,k is then obtained from equation 



([9]) by plugging-in equations ([2]) and ([3]), where e(fk,r) after 
decoupling is numerically obtained at any frequency from 
equation ([H]). 



2.3 Pulsar timing observability 

We are interested here in MBHBs potentially observable 
with ongoing and forthcoming PTAs. At the leading or- 
der, an eccentric MBHB radiates GWs in the whole spec- 
trum of harmonics f r , n — nf r ,k (n = 1,2,...). The sky- 
polarization averaged amplitude observed at a frequency 
f n = ,fr.n/(l + z) a t comoving distance d is given by 
|Finn fc Thorne!l2000r i 



hn(fn) 



32 M 5/3 
5 nd 



(2nf r , k ) 2/3 ^(^e) 



(10) 



where M = M 1 3/5 M 2 3/5 /(Mi + M 2 ) 1/5 is the chirp mass of 
the system and g(n, e) is a combination of Bessel functions 
that quant ifies the relative power radi ated in each single har- 
monic fsee lAmaro-Seoane et alj|20ld . for a detailed descrip- 
tion). If we observe for a time T b s , the relevant detectable 
amplitude of each harmonic is approximately 



'iobs.n — % V lobs/ii 



(11) 



where \/T ^ s f n is simply the number of cycles completed by 
the n-th harmonic in the observation time. The timing resid- 
uals are defined as inte grals of the GW dur ing the observa- 
tion time. As detailed in lSesana et al.l i|2009l ). each harmonic 
induces an average timing residual of the order: 



^Cgw\/nJ — 



8 /lob 



15 2nf n 



(12) 



where -y/8/15 is the average 'antenna beam pattern', i.e. 
the average of the signal performed over all possible source- 
pulsar relative orientations. The total residual can then be 
assumed to be of the order: 

1/2 



St, 



/ j "t gw (/n / 



(13) 



In observations with PTAs, radio-pul sars are moni tored 
weekly for total periods of several years i|Hobbsll2011f) . As- 
suming a repeated observation in uniform At time intervals 
for a total time Tobs, the maximum and minimum resolvable 
frequencies are / max = l/(2At) m 10 _6 Hz, corresponding to 
the Nyquist frequency, and / m i n = 1/T bs ~ 3 x 10~ 9 Hz for 
a 10 yr observation time. We are therefore interested only 
in MBHBs producing a significant residual in this frequency 
window. 

Coupling the cosmological coalescence rate to the de- 
tailed binary evolution (see Sections 2.1 and 2.2) we numer- 
ically obtain a distribution d 4 N/ (dMidM2dzdlnf ri k) so that 
J d 4 N / (dMidM2dzd\nf rt k) in the appropriate mass, redshift 
and frequency ranges, give the 'average' number of sources 
observable in the Universe taking an ideal snapshot of the 
whole sky at any time. To quantify the population of relevant 
PTA-MBHBs, we generate Montecarlo samples of MBHBs 
according to such distribution assuming f T> k > 10~ 10 Hz, 
M1&M2 > 10 7 M Q and z < 3. For each MBHB we then 
compute the average timing residual in the observed fre- 
quency range [/mm, /max] according to equation (fl3|) . and 
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we keep only those binaries producing a residual larger than 
0.1ns. 



3 DESCRIPTION OF THE PTA-MBHB 
POPULATION 

In our default model, we assume that MBHBs evolve ac- 
cording to the scheme described in Section \2. 21 We consider 
/3-discs, with viscosity proportional to gas pressure only and 
a vi scosity parameter (jShakura fc Sunvaevlll9731 ) a — 0.3 
(see King e t al. 20 07J, for a d etailed discussion), rh — 0.3 



ness in driving the p eriodic streaming activity. Nonetheless, 
iTanaka et al.l ( 20111 ) showed that also in this detached phase 



Kollmeier et al . 2006; L abita et al . 2009) and ef- 



(see, e.g., 

ficiency e = 0.1 (appropriate for not-to-mildly rotating black 
holes). All the relevant features of the resulting PTA-MBHB 
population, averaged over 100 Montecarlo realisations, are 
shown in figure [T] Typical PTA sources are very massive 
(M > 10 8 ), cosmologically nearby (z < 2) binaries. Mass 
ratios are in the range 0.1-1, with a long tail extending 
to 10 -3 . Given the significant eccentricities, orbital peri- 
ods contributing to the PTA signals extend up to > 100 
years. For the range of masses involved, this corresponds to 
a broad distribution in orbital semimajor axis in the range 
30 — 10 3 -Rs and circular velocities peaked around 10 4 km 
s _1 (if the systems were in circular orbits). Typical coa- 
lescence times (t c , defined as the integral of equation ((9} 
from the restframe binary frequency to coalescence) are be- 
tween 10 4 - 10 6 years (bottom left panel); PTA-MBHBs are 
therefore caught in the final few hundred thousand years of 
their life. The bottom-central panel shows the distribution 
of induced residuals. Order of <; 1000 sources contribute to 
a level of Ins or more, which can realistically be considered 
the ultimate goal of future PTA efforts. Their superposi- 
tion will result in a confusion-noise-type foreground, similar 
to the WD- WD binar y signal expected for LISA (see, e.g., 
iNelemans et al.ll200ll ). As in this latter case, few (maybe 
most) of the sources contributing to the signal may be in- 
dividually resolvab le. Preliminary, pu rely frequency based, 
crude estimations l|Sesana et al.ll2009r ) forecast resolvability 
of at least ~ 5 — 10 sources, but many additional sources 
will likely be resolved using the spatial information enclosed 
in the array (up to 2 N/7 per frequency bin , for an array 
of N pulsars according to lBovle fc Penll2010l . the prescrip- 
tion we will use below). The bottom-right panel shows a 
rough estimation of the mean X-ray flux on Earth assum- 
ing that the accretion rate of the binary (fueled by streams 
flowing through the gap, see next section) is the same as 
the one assumed for modeling the outer circumbinary disc 
(i.e., rh = 0.3 in this case) and a bolometric correction of 
3% (|Lussoll2010l ). Given their high mass and low redshifts, 
PTA-MBHBs are quite bright, with fluxes generally above 
10 13 erg s" 1 cm -2 . In all panels, blue-solid histograms re- 
fer to sources at a > a fl . In such systems, the gas can fol- 
low the binary during its inspiral, and strongly periodic in- 
flows of gas feeding the t wo black holes are a robust predic- 
tion of SPH simulations (jRoedig et al.ll20Hl . see next Sec- 
tion). Thus, we can safely count them in the population of 
PTA-MBHBs suitable for multimessenger astronomy obser- 
vations. The dashed-red histograms, instead, track binaries 
at a < a fr . For such systems, torques exerted by the binary 
quadrupole moment on the inner edge of the disc become 
smaller and smaller as the binary shrinks, losing effective- 



(which we labelled phaselll) gas can efficiently leak through 
the gap, feeding the central binary; moreover, residual activ- 
ity related to the consumption of the fossil ga s left around 
the t wo MBHs can still be present (see, e.g., IChang et al.l 
[2010J). Therefore, detached systems should also be interest- 
ing targets for electromagnetic counterpart identification; 
however, we do not consider the m in the following d iscus- 
sion, and we refer the reader to ITanaka et al.l (|201lh for a 
comprehensive analysis of the subject. 



3.1 Parameter dependence and caveats 

The populatio n shown in figu re [T] relics on a selected MBH- 
host relation (IGiiltekin et alj 120091 ) and on a specific disc 
model (the /3-disc model) with fixed parameters. Even 
though the qualitative features of the predicted popula- 
tion are robust against our particular choices, the effec- 
tive number of sources and the relative population of cou- 
pled/decoupled systems are not. Both lowering a and rh re- 
sult in longer migration timescales in the disc. This, on one 
hand, increases (up to a factor of a few) the number of PTA- 
MBHBs. However, in this case, GW emission takes over at 
slightly larger separations (see the weak dependence of a r on 
a and rh in equation (O), and more binaries (especially with 
periods < 5 years) will be already detached from their discs 
and probably less suitable for multimessenger observations. 
In an a-disc, the viscosity is much larger in the radiation 
dominated zone, and, even though there is no self consistent 
solution to the binary-disc evolution, it is likely that the de- 
coupling would happen at smaller separations. In any case, 
t v is surely much smaller in the radiation dominated zone 
for an a-disc, and the disc 'freezing' radius would certainly 
be smaller, allowing binaries to be in touch with their cir- 
cumbinary discs for longer time and at smaller separation. 
In this respect, the ,5-disc assumption used here is likely to 
be conservative. The impact of the chosen MBH-host rela- 
tion is also mild, affecting the PTA-MBHB population by 
a factor of two-three at most. Also the binary eccentricity 
is not very important in terms of the population itself. In 
our simple model, in fact, we assume the disc migration to 
be independent on the binary eccentricity. A population of 
circular binaries will only decouple at mildly smaller val- 
ues of a, slightly increasing the population of coupled sys- 
tems. However, in this case, variability related to periodic 
streams flowing thr ough the gap would be highly reduced 
IjRoedig et alj|201ll ). making probably difficult to recognize 
these systems through periodicity studies. 

As a note of caution, we remind the reader that our 
models assume that all merging MBHBs are active in their 
last evolutionary stage. Under such condition (and the as- 
sumption that systems are accreting at rh — 0.3), the source 
flux distribution shown in figure \T\ accounts for ~ 30% of 
the s oft-X AGN LogA-LogS as observed by ROSAT l|Voges! 
Il999i ). meaning that one bright X-ray AGN out of three 
is indeed a massive binary. This, although not inconsis- 
tent with any observation, is, of course, a strong statement, 
that can be relaxed if rh < 0.3, if a certain fraction of 
our sources is indeed obscured, or simply if only a frac- 
tion T of MBHBs in indeed active. In this latter case, it 
is sufficient to rescale our predictions by a factor _F, ac- 
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Figure 1. General features of the typical MBHB population contributing at a level of 0.1ns or more to the PTA signal in the 3 X 10 — 
10~ 6 Hz frequency window. In each panel the blue-solid histogram refers to sources at a > a , the read-dashed histogram refers to sources 
at a < a' r (see text for details). In the top row, we plot the cosmological evolution-related MBHB properties: the primary mass M%, mass 
ratio q and redshift distributions, from the left to the right. In the central row we plot MBHB properties describing the dynamics of the 
system: the binary semimajor axis in units of Rg, the binary period and the circular velocity V c distributions, from the left to the right. 
In the bottom row we plot the coalescence time distribution (left), the distribution of the induced timing residuals R (center) and the 
X-fiux (0.5-10 keV) assuming a bolometric correction of 3%. Distribution refers to the average over 100 Montecarlo realisations of our 
default model: a = 0.3, m = 0.3, e = 0.1, and e^ec given by equation (71 . 



counting for the fraction of active binaries. Similarly, we as- 
sumed that all MBHBs overco me the 'last parsec problem' 
IjMilosavlievic fc Merrittll200lT ); which, in the light of several 
recent results about both gas and stellar driven MBHB hard- 
ening (|Escala_jrt al. 2005: D otti et al.; 2007; Cuadr a et al.l 
l200a : IColpi et alj|2009l ; iKhan et al.ll201ll : |Preto et al.ll201lf ). 



we consider a reasonable assumption. We will see in the 
following that even for T < 0.1, the number of predicted 
sources is sizable, making multimessenger astronomy with 
PTA sources an appealing prospect, worth of deeper inves- 
tigations. 
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4 DYNAMICAL MODELLING OF THE 
PROTOTYPICAL PTA-MBHB 

Having studied the MBHB population relevant to our inves- 
tigation, we model in this section the dynamics of the pro- 
totypical PTA-MBHB surrounded by a circumbinary thin 
accretion disc. We use a syncretic approach, combining an- 
alytical accretion disc models to high resolution SPH simu- 
lations. The former provide the global set-up of the system, 
the latter are necessary to gather a better insight of the dy- 
namics of the torqued material streaming through the inner 
edge of the circumbinary disc. 

4.1 Standard model for a PTA binary 

We model a MBHB of a total mass M = Mi + M 2 sur- 
rounded by a thin Keplerian gaseous disc of mass Ma- Al- 
though Newtonian simulations are in principle scale-free, the 
expected properties of the circumbinary disc depend on the 
assumed mass and semimajor axis of the binary. We there- 
fore need to work out an a priori scaling of the SPH simu- 
lation. We set primary mass M\ = 2.6 x 1O 8 M0, mass ratio 
q — M2/M1 — 0.35 and semimajor axis ao = 0.012 pc, con- 
sistent with the typical range of source param eters found in 
the pr evious section (see figure [l}. Following iRoedig et al.l 
(|2011f ). we set the the initial binary eccentricity to a value 
eo = 0.6. The binary has angular velocity flo = (GM/ajj) 1 ' 2 
resulting in a period P = f^ 1 — 27r/f2o w 6 years (/o is the 
orbital frequency). The binary is surrounded by a corotat- 
ing coplanar thin disc, a configuration expected in the late 
stage of a relatively ' wet' galaxy merger, as shown by sev - 
eral SPH simulations JMaver et alj|2007l : iDotti et al.lf2 009). 
A tru ncated circumbinary /3-disc has mass IjHaiman et al.I 
I2009D 



M d = 8.63xlO 4 M aj;3 /5 



■m . :i 
eo.i 



7/10 



M 8 11/5 (i?; 



5/4 _ 
out 



< /4 ),(14) 



where two limiting radii of the disc, R[ n and R ou t, are ex- 
pressed in units of 10 3 -Rs and, in our simulation, correspond 
to Ri n = 2ao and i? ou t = IO00, respectively. Plugging in our 
default disc model (ao.3 = rko.3 = eo.i = 1) an d the binary 
parameters assumed above, we obtain R[ n w 0.6, -Rout ~ 3 
and M d » 5 X 10 6 M Q « 1.5 x 10" 2 M. Comparing the vis- 
cous timescale of the binary to the GW driven orbital de- 
cay (equation JTJ and @, in Section 2.2), we find that the 
MBHB is still attached to the circumbinary disc, and we 
can study the dynamics in the Newtonian approximation, 
neglecting relativistic effects. 



4.2 Numerical realisation 



We use the SPH cod e Gadget-2 l|Springell 12005) as done in 
IRoedig et alj (|201ll ). we refer the reader to that paper for 
details and only outline the essentials here. As initial con- 
dition, we use the binary outlined in the previous section 
surrounded by an 8 million particle circumbinary disc, and 
let it relax over 9 orbits. The dynamics of the MBHs, which 
are modelled as Newtonian point masses, is followed with 
a fixed time-step, equal to 0.01 fig" . The gas is allowed to 
cool on a time scale which is proportional to the local dy- 
namical time of the disc thus we set f3 — t coo i/td y n = 10. 
In order to confine the gas in the central cavity to a thin 



geometry, we assume that the small amount of gas present 
in the inner cavity (r<1.75a) is isothermal, with an inter- 
nal energy per unit mass u « 0.14(GAf/_R). The softening 
for the particles is adaptive with a minimum of 0.001 ao, 
while the MBHs are not softened Q The sink radius, the 
radius b elow which any passing particles is counted "ac- 
creted" IjBate et alj|l995f ). is set to 0.005 ao, which corre- 
sponds to w 6 Rs for the secondary MBH. The necessity of 
such a high resolution is given by the high mass ratio be- 
tween holes and discs. Furthermore, we need to be able to 
resolve sub-Eddington accretion rates onto the black holes, 
which is MBdd = 8M0/yr, thus we need to have single par- 
ticles that are some factor smaller than this. In our choice of 
units, the mass of a particle is Msph — 0.79 Mq . The aspect 
ratio of the disc is h/r ~ 0.04 with a mass ratio between disc 
and holes M d = 1.58 x 10~ 2 M. 



4.3 Variable inflows and minidisc formation 

Gravitational torques, acting on the inner edge of the cir- 
cumbinary disc, cause periodic inflows of mass that are trig- 
gered by the secondary pulling gas out of the disc edge at 
each apo-apsis passage. We follow the streams down to the 
sink radii of the two black holes and monitor the particles 
that are swallowed by the MBHs and bin them according 
to the timestep when they were swallowed. An example of 
the forming streaming structure is given in figure [5] where 
we show the linear column-density of the gas in the gap 
Q The inflows are well resolved, and the gas bound to each 
MBH is visible to form a sense of minidiscs. Their radial size 
is about 0.05 ao and their average density p — 2 x 10~ 5 g 
cm -3 , which gives a total mass of the minidisc m ln d 1 ~ 3 M© 
for the primary and m m d 2 ~ 0.5 Mq for the secondary. The 
detailed structure of the inflows and their physical proper- 
ties strongly depend on the disc thermodynamics; moreover, 
their behaviour close to the MBHs is likely to be affected by 
the comparable size of the sink radius, the smoothing length 
and the physical size of the minidiscs. Nonetheless, we can 
consider the rate at which particles cross the sink radii of 
the two black holes (what we call numerical accretion rate) 
as a good proxy to the rate at which material is fed from the 
outside circumbinary disc to the inner region. Such a rate 
is shown in figure [3] over sixty orbital periods. The orbital 
periodicity is quite striking, with a maximum-to-minimum 
ratio larger than a factor of three. The maximum accretion 
rate declines from m ~ 1 to m ~ 0.4, as a result of the relax- 
ation of the disc from its initial condition which is slightly 
out of equilibrium. We notice that the inflows occur on the 
MBHB orbital period, and are triggered by the secondary 
pulling material off the inner edge of the circumbinary disc 
at its apo-apsis passages. Their nature is therefore primar- 
ily gravitational, and should be weakly dependent on the 
exact thermodynamics of the disc. Our simulation therefore 
suggests that gravitational torques are effective in feeding 
periodically the central gap at a rate which is a substantial 



3 Note that this does not introduce any artificial effect in the com- 
putation of the gravitational forces, since the MBHs are modelled 
as sink particles with sink radius larger than the typical softening 
of the gas particles. 

4 fig.[2]made using SPLASH l|Pricell2007f) . 
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og column density 

Figure 2. Logarithmic column density of the cavity, the MBHs 
are inserted as small black dots, the axes are in units of the semi- 
major axis an 



fraction of the Eddington rate. Our results are in agreement 
with similar full GR hydro simulatio ns studying the ev ofu- 
tion of MBHBs at closer separations (|Farris et al.ll201ll ). in 
which similar streams feeding the two MBHs are observed 
down to « 10M (r^3x 10" 4 pc for our MBHs). 



4.4 Analytical modelling of the minidiscs 

As stated in the previous section, the detailed behaviour of 
the inflows close to the MBHs is likely to be affected by the 
comparable size of the sink radius, the smoothing length 
and the physical size of the minidiscs. We therefore decide, 
for the purpose of the following discussion, to build a sim- 
ple analytic parametric model for the minidiscs. We assume 
that the infalling material is fed from the outer disc at some 
given rate rh. To simplify the notation, we do not distin- 
guish between the two MBHs; all the quantities appearing 
in the equations below, are intended to be specific of the 
MBH under consideration^, and the results apply to both 
MBHs (and both attached minidiscs). The infalling mate- 
rial is captured in the Roche lobe of the i-th MBH with an 
average specific angular momentum J with respect to that 
MBH. This is converted to a characteristic size of a Kep- 
lerian orbit around the MBH simply as r m d = J 2 /GM. If 
r m a < 3Rs the inflow is radial and accretion proceeds in a 
Bondi-like fashion. Otherwise, the infalling material settles 
at a characteristic radius r m( j around the MBH and dis- 
sipates its angular momentum through viscous processes, 



5 For example, M s = M;/1O 8 M is the mass of the i-th MBH in 
units of 10 8 solar masses, an. 3 = Oi/0.3 is the viscosity parameter 
of the i-th disc in units of 0.3, and so on. 



Figure 3. Numerical accretion rate as a function of time (see text 
for details). The periodicity is remarkable, as also highlighted in 
the upper-right zoom-in. 



diffusing inwards and eventually being accreted by the two 
holes (a process similar t o the debri s fallback in tidal dis- 
ruption induced accretion. lReeslll98a ). The time needed for 
an annulus of material at radius r to be accreted is 



M & r[ 



7/2 



-2/5 



M 



6/5 7/5 



(15) 



(16) 



. n Q1 -r (mo.3 

iaco,a = 0.31 yr a 3 

\ eo.i 

i coo - 4 /5 I 1TI-0.3 

tacc,/3 = 683yra 3 ' 

V eo.i 

for a and /3 discs respective!^. Here r\ is the size of the 
minidisc in units of WRs- If iacc is longer then the binary 
orbital period P, then a pair of persistent, periodically fed 
minidiscs is formed around the two MBH. If, conversely, i acc 
is shorter then P, then the streaming periodicity is reflected 
in episodes of periodic accretion. For o-discs, this happens 
for 

4/7 

M 8 - 2/7 P y 2 / 7 (17) 



. on D 2/7 frho.3 
r m d,max < 20 its a .3 ' — 



V £ o.: 



(P yr is the MBHB period in years), while it never happens 
for /3-discs. Figure U shows the maximum radius r m<1 ,max 
and the corresponding maximum mass m m d,max of an a- 
minidisc. For average feeding rates of rh w 0.3 (consistent 
with the output of our SPH simulation) , maximum disc sizes 
are in the range 10-40 Rs- 

In our specific simulation, the average angular momen- 
tum of the particles captured by the MBHs within their re- 
spective Roche-lobes implies r m d ~ 10— 20Rs, meaning that 
the accretion process must happen through viscous inspiral 

6 To sketch the situation, we describe the minidiscs as small 
steady accretion discs, even though the condition of stationar- 
ity of the system are not met because of the periodicity of the 
infalling streaming. 
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Figure 4. Characteristic properties of Q-minidiscs. Upper left 
panel: maximum minidisc radius r m( j max allowing for t a cc,a < P, 
as a function of P. Thick-red curves are for m = 0.3 and a = 0.3 
and thin- blue curves are for m = 0.1 a = 0.1. Short-dashed, 
long-dashed and solid curves are for a MBH mass of 10 8 , 10 85 
and 10 9 Mq respectively. Upper left panel: maximum mass of the 
minidiscs enclosed in r m d, m ax as a function of P, line and color 
style as in the upper left panel. Lower panel: maximum size of 
the minidiscs relative to a as a function of P, compared to the 
Roche- lobes of the two MBHs (horizontal ticks). Here we assume 
our default binary; long-dashed lines are for Mi, short-dashed 
are for M2. 



rather than through radial infall, disfavoring the Bondi-like 
accretion scenario. Moreover, r m( j < r m d, ma x for both MBHs; 
for Q-discs, this suggests MBHB period-related periodicity 
in the accretion and thus in the emitted luminosity. 



5 EMITTED SPECTRUM AND 

ELECTROMAGNETIC SIGNATURES 

In optically thick disc theory, each disc annulus emits black- 
body radiation corresponding to its particular temperature; 
for a standard opt ically thick geometrically thin disc, this is 
JFrank et alj|2002h 



T(r) = 



3GMMf(r) 



87rr 3 cr s b 



1/4 



(18) 



where f(r) — (1 — SRs/r) 1 ' 2 (we assume a Schwarzschild 
black hole), and <7 s b is the Stephan-Boltzmann constant. The 
emitted luminosity at frequency v is given by the integral 
over the disc extension of the blackbody radiation emitted 
by each annulus, i.e. 

32nhpu 3 r° ut rdr 



L v = 



ch P u/kT(r) 



(19) 



where 7m and r ou t are the inner and outer edges of the 
disc respectively, hp is the Planck constant, and k is the 
Boltzmann constant. In our model, each of the discs (the 



circumbinary plus the two minidiscs) emits a thermal com- 
ponent according to equation (|19fl . We therefore infer the 
presence of three distinctive continuum components: (i) a 
thermal component peaked in the optical/IR emitted by 
the circumbinary discs; (ii) two thermal component peaked 
in the UV associated to the inner minidiscs; (iii) two X- 
ray powe rlaw associated to t he tenuous hot electron plasma 
(corona, iGaleev et al.l 1 19791 ) surrounding the inner mini- 
discs. These latter components are generated by the in- 
verse Compton scattering of UV photons emitted by the 
inner radii of the minidiscs against the diffuse hot plasma 
of electrons embedding the inner part of the discs. The up- 
scattered photon spectrum can be modelled as a power-law 
spanning the soft/hard-X domain. For the qualitative na- 
ture of our di scussion, it is suffi cient to consider a flat spec- 
trum in vL v (JHaardt fc M araschi 1993), normalized to give 
£o.5-iokcV = 0.03Lboi (|Lussoll2010l ). The general features 
of the predicted continuum are shown in figure [S] for two 
selected systems. The presence of a gap is reflected by the 
characteristic double bump; the one in the IR is due to the 
circumbinary disc, while the one in the UV is produced by 
the two inner minidiscs. 

On top of the continuum, several sets of emission lines 
are also expected, as in common AGNs. In particular: (i) 
optical broad emission lines caused by the photoionization 
of the tenuous infalling material and the outer circumbi- 
nary disc (r > O.Olpc) by the inner ionizing UV source (the 
minidiscs, r < lCF 3 pc); (ii) two 6.4keV fluorescence iron 
broad emission lines (Fe Ka lines) produced by the reflection 
of the up-scattered corona X-ray photons on the surface of 
the inner accretio n minidiscs at only few Schwarzschild radii 
IjMatt et al.lll991f ). Moreover, several other spectral features 
may be present: (i) some weak op tical/UV emission asso- 
ciated to the str eaming material (|Bogdanovic et al.l 120091 ; 
iDotti et al.ll2009l ); (ii) X-ray hot spots related to the in- 
streaming material shocking onto the outer edge of the inner 
minidiscs (if t acc > P, see Section 5.1). 



5.1 Characteristic signatures 

Depending on the strength and angular momentum of the 
inflows, the features outlined above may or may not be 
present. We identify three different scenario: 

(i) 7"md < 3Rs- The streams flow radially onto the two 
MBHs. The radiation emission of the associated Bondi- 
like accretion is likely to be extremely inefficient. No ion- 
izing UV continuum and, consequently, no optical broad 
lines are present. In this scenario, no signatures of the 
existence of an accreting system may be detectable. The 
only guaranteed component is the continuum associated to 
the circumbinary disc thermal emission. If, however, the 
spread in specific angular momenta of the accreted mat- 
ter is large compared to the angular momentum of the last 
stable orbit (Aj T > Rsc), then a scenar i o ana log to the 
one proposed by llllarionov fc Beloborodovl (120011 ) for wind- 
fed high mass X-ray binaries may take place. The inflowing 
streams collide in a caustic at few Rs and are promptly ac- 
creted in a free-fall time, forming a small scale, inviscid, ra- 
diatively efficient accretion disc i|Beloborodov fc Illarionovl 
l200ll ; IZalamea fc Beloborodovll20091 ). In such models, the ac- 
creted matter/luminosity conversion efficiency ranges from 
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Figure 5. Continuum components of the emitted spectra. In each 
panel we plot the thermal emission of the circumbinary disc (blue 
dot-dashed line), and of the two inner minidiscs (long-dashed red 
line for the primary, short-dashed green line for the secondary); 
depending on the nature of the inflowing streams, these latter 
components may or may not be present. The black solid line is 
the resulting total continuum, where a flat X-ray emission from 
a putative hot corona has been added (see text). The vertical 
dashed line is the energy center of a typical green optical filter. 
The two panels refer to two fiducial systems with parameters 
labelled in figure. In the upper panel, each minidisc has size r md = 
25R$ of the correspondent MBH; while in the lower panel, the 
size of the minidiscs is r md = 15-Rg only. We assumed our default 
accretion model (a = 0.3, rh = 0.3, e = 0.1). 



0.03 for Schwarzschild MBHs, up to 0.1 for maximally spin- 
ning MBHs. The MBHB would therefore appear as a peri- 
odically variable, luminous X-ray source. 

(ii) 3Rs < r m( j < -Rmax, Q-discs only. The streams form 
standard minidiscs around the two MBHs. Their typical vis- 
cous time is shorter than the binary orbital period; therefore, 
accretion onto the two MBHs can be directly linked to the 
periodic streams. In this scenario, efficient thermal emission 
comes from both the circumbinary disc and the two mini- 
discs. While emission from the circumbinary disc will be 
relatively steady, strong periodicity in the UV/soft-X will 
reflect the periodic accretion related to the minidiscs. The 
resulting periodic ionizing photon flux, will cause optical 
broad line strength to fluctuate by a factor larger than two 
over the orbital period. However, identification of such opti- 
cal variability may be challenging, since it would require an 
all sky periodic spectros copic monitoring. Large photometric 
surve ys like the LSST IjLSST Science Collaborations et al.l 
12009 ). will be sensitive to the luminosity of the optical con- 
tinuum, which, depending on the characteristics of the sys- 
tem, may be dominated by the steady circumbinary disc (see 
figure [5]). Nonetheless, depending on the amount of redden- 
ing, color selection may identify candidates showing a dis- 
tinctive optical bump across the different wavelength bands. 
Note that in such close systems, the broad line emission re- 



gion is likely related to the circumbinary disc; the typical size 
of the broad line emission region in AGN is ~ 0.01 — O.lpc, 
much larger than the size of the minidiscs, and correlates 
with the AGN luminosity JKaspi et al.ll2005l ). Detection of 
separate sets of broad optical emission lines related to accre- 
tion onto the individual MBHs, or the detection of a broad 
emission line shifted with respect to the narrow em ission 
line associated to the host (see, e.g. iDecarli et al.ll2010l ) are 
therefore unlikely in this case. Periodic soft-to-hard X-ray 
emission in response to the periodic continuum may be as- 
sociated to the hot corona, and variable relativistic Fe Kq 
lines may also be present. This is the scenario suggested by 
our SPH simulation. 

(iii) r md > -Rmax (a-discs), or r md > 3Rs (/3-discs). 
The streams cannot be swallowed by the MBHs within an 
orbital period. Steady thin accretion minidiscs form around 
the two MBHs. The accretion rate may be still quite vari- 
able, but redistribution of angular momentum within the 
minidiscs makes impossible to link the accretion variabil- 
ity to the periodicity of the streams. Two superposed rel- 
ativistic Fe Ka lines associated to the minidiscs are likely 
to be observable. Moreover, shocks created by the collision 
of the instreaming material and the minidiscs at their outer 
edge may result in X-ray hotspots varying on the orbital 
period. The luminosity and electromagnetic frequency at 
which these regions irradiate depend on the relative velocity 
(Urei) of the streams with respect to the minidiscs. We can 
estimate v IC \ for M2, when the secondary is at the apocen- 
ter, assuming that its minidisc fills the whole Roche lobe. 
In this case, the emission is expected to peak in the hard 
X-ray, at about ~ v^ cl fi p /k>10 keV (fi p is the mean molec- 
ular weight of the plasma). The luminosity depends on the 
time-scale r s hock over which the shock irradiates its internal 
energy. Assuming r s hock to be of the order of the orbital pe- 
riod of the gas at the edge of the minidisc (Pdisc), we obtain 
ishock « mv^ cl P/ P disc ^10~ 3 L Edd . Smaller r shock and larger 
i>rd result in larger luminosities, but increasing v rc i shifts the 
peak of the emitted luminosity at considerably higher fre- 
quencies. Note that low values of v IC \ are the most favorable 
conditions to efficiently perturb and bind large amounts of 
gas. 

In the following we will focus on two distinctive signa- 
tures in the X-ray domain. 

• Periodic X-ray emission related to the periodicity of the 
minidisc fed through the gap. The outer circumbinary disc 
is relatively stable, and optical variability related to the 
streams may be overwhelmed by the outer disc continuum. 
Periodic variability, related to (i) periodic accretion of the 
minidiscs, (ii) the associated varying flux of UV photons up- 
scattered in the corona, (iii) X-ray hotspots created by the 
periodic streams colliding with the minidiscs, may instead 
be easily detectable in UV and in X-ray, 

• Double relativistic fluorescence Fe Ka lines. Such 
lines are a common feature i n AGNs IjNandra et al.l 120071 ; 
Ide La Calle Perez et al.l 120101 ) , and are produced at only 
few Schwarzschild radii. The line profile strongly depends 
on the location of the last stable orbit orbit around the 
MBH (which depends on the spin magnitude) and on the 
disc inclination with respect to the observer. If accretion is 
efficient on two MBHs with different spin parameters (mag- 
nitude and/or orientation), a distinctive 'double Fe Ka line' 
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feature may be observable in the hard X spectrum of the 
source. 

In the two following sections we will consider these 
two possibilities separately, proposing possible strategies for 
combining X-ray and PTA observations in the final discus- 
sion. 



6 X-RAY PERIODIC VARIABILITY 

6.1 Sub-population of observable periodic sources 

Being interested in variability related to the binary orbits, 
only a subset of the population shown in figure[T] having rea- 
sonably short periods will be suitable targets. We therefore 
limit our study to the systems with P < 5yr. The result- 
ing population, represented in the left plot of figure [6] for 
our default model, retains the mass, mass-ratio and redshift 
distribution of the parent one. Detached systems obviously 
have a period distributions extending to P < lyr, while 
binaries still attached to their discs (relevant to our investi- 
gation) are abundant at P > lyr. Assuming that a fraction 
of 3% of the bolometric luminosity is emitted in the X-ray[j, 
we can construct the LogiV-LogS' (i.e., the cumulative num- 
ber of sources having a measured 0.5 — lOkeV flux larger 
than a certain value S, as a function of S) of the emitting 
population. This is shown in the right plot of figure [6] for 
both our default population (upper panel), and for an alter- 
native model in which the smaller a viscosity and accretion 
rate imply an earlier detachment of the circumbinary disc. 
At large luminosities, such periodic binaries may add up to 
5% of the observed luminous X-ray sources in th e sky. The 
curre ntly operating X-ray all sky monitor MAXI (jMatsuokal 
120091 ) may already have detected an handful of them in its 
surveys. On the other hand, according to our default model, 
there might be up to 500 MBH binaries significantly con- 
tributing to the PTA signal, showing periodic variability on 
a timescale of 1 — 5 years, in the sensitivity reach of the up- 
coming eRosita observatory i|Predehlll2010r ). We stress, once 
again, that our models assume that all merging MBH bina- 
ries are accreting in their last evolutionary phase prior to 
coalescence, which may be in fact likely in gas rich merg- 
ers, but nonetheless it remains an ad hoc assumption. If, 
for instance, only 10% of the merging systems is active, the 
number of sources showing emission periodicity drops to 10- 
50, which is however still an interesting, sizable sample. 



6.2 Simulating observations: sampling and 

statistics 

An X-ray luminosity above the instrument sensitivity is 
certainly not enough to identify a periodic source. The 



7 This picture is appropriate for hot corona reprocessing. Inviscid 
minidiscs have matter/luminosity efficiency conversions similar 
to standard accretion discs (e ~ 0.03 — 0.1), and even though 
detailed calculation of the spectrum are not available, we may 
expect comparable X-ray luminosities. If the minidiscs are instead 
steady, then X-ray luminosity emitted by shock induced hot spots 
may be more than an order of magnitude smaller (10 — 3 LEddi 
compared to 3% of 0.3L.Edd f° r corona reprocessing in our default 
model). 



lightcurve has to be reconstructed with a minimum num- 
ber of datapoints per orbit, for at least few orbits. In this 
subsection we carry a detailed statistical study of the period- 
icity observability. The main goal is to provide the minimum 
required sampling cadence an observatory must have to effi- 
ciently identify periodicities. Since astronomical time-series 
are always limited to a certain amount of data points, we 
address the problem of having a very coarse sampling rate of 
our lightcurve. We assume that the X-ray lightcurve mimics 
the numerical accretion rate of figure [3] 

The complete lightcurve of our simulation consists of 
L — 898 points, sampling ~ 80 orbits; each orbit is therefore 
sampled by ~ 11 equally spaced points. We create subsam- 
ples of this lightcurve by selecting a random starting zero 
point p and then considering all the subsequent np points, 
where n — 1, ..., N — 1, and p £ [1, 2, 3]. A stride of p = 1 
translates into taking 11 points per orbit (points p, p + 1, 
p + 2...), p — 2 means taking ~ 6 points per orbit (points 
p, p + 2, p + 4...) and p = 3 means taking ~ 4 points per 
orbit (points p, p + 3, p+ 6...). We refer as GPopl to a sam- 
ple covering six orbits with all points included in the sam- 
pling (total of N — 65 points); 3Pop2 refers to a string of 
three orbits sampled every other point (total of N = 18 
points), etc. The reason to consider different sampling is 
to investigate the effect of the observation cadence. Assum- 
ing a binary with, e.g. 3 year period, pi corresponds to a 
three month cadence, p2 to six months, etc. If the number 
of points covered by the typical subsample is N, we can draw 
a maximum of f = L — N subsamples Q To account other 
possible sour ces of variability (A GNs are generally variable 
in the X-rav. lGrupe et al.ll200fi ). the samples are then pol- 
luted with Gaussian noise H with a dispersion a = 10% 
of the maximum accretion rate/luminosity of each specific 
sample. Each sub-sample is Fourier-tra nsformed usin g the 
Normalized Lomb-Scargle Periodogram (|Scarglelll982r i and 
the False-Alarm-Probability (FAP) of the highest detected 
peak is computed. After a cross-check that this highest 
peak corresponds to the fundamental frequency of the bi- 
nary /o = 1/Po, the FAPs of the /o frequency are binned 
logarithmically and shown in the histograms in figure[7] The 
left panel shows that a minimum of ~ 15 — 20 datapoints 
covering at least three orbits is needed in order to iden- 
tify a decent fraction ~ 50% of the sources at a 10% false 
alarm probability level. Increasing the number of datapoints 
dramatically increases identification performances; already 
with ~ 30 points in the lightcurve we can detect most of the 
sources ~ 70% at a false alarm probability level of 1%. 

Additionally, it is of interest, in order to crosscheck with 
a possible PTA identification of the source, to which accu- 
racy the fundamental frequency can be recovered from the 
Fourier Spectra. For an observation time T, the frequency 
resolution bin of the observation is 1/T. If therefore we ob- 
serve a frequency /o, the relative accuracy of the observa- 
tion is A/0//0 = (l/T)//o. However, we can oversample 
the spectrum to get a better constrain on the observed fre- 



8 Note that the subsamples in general overlap with each other, 
therefore cannot be regarded as uncorrelated. However, a much 
longer simulation would be needed to draw a large number of 
uncorrelated samples, which would have been prohibitively time 
consuming. 
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Figure 6. Left plot: properties of the PTA-MBHB population observable through periodicity, contributing at a level of 0.1ns or more to 
the PTA signal in the 3 X 10 — 9 — 10 — 6 Hz frequency window in our default model, averaged over 100 Montecarlo realisations. Systems 
with orbital periods < 5 years only were selected. From the top left to the bottom right, we plot the primary mass Mi, mass ratio q, 
redshift and period distributions. Right plot: periodic MBHB contribution to the X-ray LogTV-LogS function, assuming a bolomctric 
correction of 3% (0.5-10 keV). The top panel is for our default model; the bottom panel is for a less optimal model (m = 0.1) for which 
MBHB-disc decoupling occurs earlier. Vertical dotted lines depict the flux limit for a single eRosita pointing and for the one-month sum 
pointings of MAXI. Dashed thin line depicts the upper end of the AGN LogTV-LogS as observed by ROSAT ( Vogcs 1999). In both plots, 
lincstyle as in figure [T] 
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Figure 7. Detection statistics of periodic sources. In each plot, we show histograms of the number of samples correctly identified within 
a certain False- Alarm-Probability (FAP). Left panel: statistics is constructed with 6 points per orbit over 3 orbits (3Pop2, red histogram), 
or with 3 points per orbit over 6 orbits (6Pop4, blue histogram). Central panel: statistics is constructed with 11 points per orbit over 
3 orbits (3Popl, red histogram), or with 6 points per orbit over 6 orbits (6Pop2, blue histogram). Right panel: statistics is constructed 
with 11 points per orbit over 5 orbits (5Popl, red histogram), or with 11 points per orbit over 6 orbits (6Popl, blue histogram). 



quency. For all periodograms shown, we have used an over- 
sampling factor (ofac) of 8. In evenly sampled data, we ex- 
pect the analysis to be only mildly dependent on this pa- 
rameter: For some choices of low ofac (< 4) in very short 
data sets, we observed the appearance of aliasing especially 
if high freque ncy componen ts were not suppressed (hifac > 1 
as defined in IjScarglell 19821 )). We thus concluded, especially 
in the prospect of real unevenly sampled astrophysical data, 
that the high frequency noise should be suppressed using 
very low hifac (< 1) and at the same time using high over- 
sampling (ofac > 4). For the identification of any frequency 
we set the requirement j /detected - /true] //true < 5%. In fig- 
ure [S] we show the dependence of the number of correctly 
identified frequencies on the number of points available in 



the lightcurve, for different samplings (i.e., for different num- 
ber of points per orbital period). We consider only frequency 
peaks identified at a 3 a significance or more. Note the conti- 
nuity of the p2, p3 and p4 curves (six, four, and three points 
per orbit, respectively), meaning that, as long as we sam- 
ple three or more orbits, the fraction of correctly identified 
systems depends strongly on the total number of points in 
the lightcurve and only mildly on the number of points per 
orbit (as long as this is larger than three). The discontinuity 
with respect to the pi series tells us that sampling six orbits 
with six points per orbit is much better than sampling three 
orbits doubling the observation frequency. Sampling at least 
three orbits is a minimum identification requirement, be- 
cause other sources of uncorrelated noise may wash out the 
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Figure 8. Percentage of correctly identified frequencies above 3 
a within a fractional frequency error < 5%, as a function of the 
number of points in the lightcurve. For each sampling (shown 
in different colors, see legend in figure), points from the left to 
the right correspond to three, four, five and six full orbits of the 
observed source. A 10% random Gaussian noise was added to each 
subsamplc. 



significance of the periodicity if the statistics is too poor. In 
general, having at least ~ 15 — 20 points in the lightcurve, 
sampling three or more orbits, is a minimum requirement 
for a confident identification of the periodicity. 



7 DOUBLE IRON LINES 

7.1 Relevant source population 

Spectral lines formed in relativistic accretion discs are dis- 
torted due to Doppler and relativistic effects producing a 
characteristic shape, which may be modelled to determine 
the properties of the space time ar ound the compact ob- 
ject (JFabian et al.l 11989 ; lLaonll99lT ). In particular, obser- 
vations of supermassive black holes at the center of galax- 
ies have revealed broad skewed lines, which have allowed 
constraints to be p laced on the spin of the black holes 



constraint s to be placed on the spin ol tne black holes 
dTanaka et alJll995l: iNandra et al.lfl997l . 120071 : lMillerll2007l ; 
Ide La Calle Perez et al.l 120101 ). The current status of our 
knowledge of rela tivistically broadened i ron lines from AGN 
is summarized in lGuainazzi et al.l (J201l|). 

Double relativistic Fe Ka lines are more likely to be 
present in a 'steady' environment, where material rising from 
the accretion discs (driven by heating and vertical stratifica- 
tion of the accreting material, magnetic turbulence, winds, 
or star disc interaction) has time to shape a tenuous hot 
electron plasma corona. They are therefore likely to appear 
in situations where £ acc > Po, but we do not exclude such 
possibility oth erwise. Broad Fe Ka lines appear to be com- 
mon in AGN (JNandra et al.1 120071 ; Ide La Calle Perez et"ail 
l2010l ). yet their identification requires a large number of 
collected X-ray photons, i.e. deep, targeted, time consum- 
ing observations. It is therefore reasonable to consider only 
sources that may be individually resolvable in PTA cam- 
paign s, and consequently localized in the sky to some accu- 
racy (|Sesana fe VecchiolfeoiOl ; ICorbin fe Cornisbj|2010l) . 

We therefore estimate the population of sources suit- 
able for double Fe Ka line detection by considering individ- 
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Figure 9. Properties of the individually resolvable PTA-MBHB 
population, assuming ten resolvable sources per frequency bin, 
and ten years of PTA observations, for our default model. All re- 
solvable sources contributing at a level of Ins or more to the PTA 
signal in the 3 X 10 -9 — 10 _6 Hz frequency window are considered. 
From the top left to the bottom right, we plot the primary mass 
Mi, mass ratio q, redshift and X-ray flux distributions. The ver- 
tical lines indicate the properties of the system studied in SI7.2I 
The total number of sources integrated over the blue histograms 
is ~ 20. Linestyle as in figure [T] 



ually resolvable sources only. Unfortunately, the concept of 
resolvability has not been deeply investigated in the PTA ob- 
servation context (and certainly not for eccentric binaries). 
In the circular binary case, a rough 'on e bin' rule estimate , 
provides ~ 10 bright resolvable binaries (jSesana et al.ll2009M . 
However, such estimate does not take into account the spa- 
tial information enclosed in the detection with an array of 
pulsars; two sources at different sky locations contribute dif- 
ferently in each pulsar, and their signals may be disentangled 
even i f their frequencies fall in the same bin. iBovle fe Pen! 
(J2010I ) estimated that, exploiting the spatial information en- 
closed in the signal, 2N/7 sources per frequency bin would 
be resolvable by an array of A^ pulsars. For our estimate we 
therefore (somewhat arbitrarily) pick the ten strongest GW 
sources per frequency bin, and impose a further cut at Ins 
(dimmer sources would not be individually detectable any- 
way). This leaves us with ~ 100 sources, the precise number 
being vastly independent on the details of the global popula- 
tion model, but only on the PTA observation time (assumed 
to be ten years). Such population is shown in figure [9] for 
our default model. Blue histograms represent MBHBs still 
attached to their circumbinary discs (i.e., with a > a r ), and 
therefore plausibly hosting small accretion minidiscs. We are 
left with ~ 20 bright, low redshift sources. 

Notice however, that our estimation of a r is quite con- 
servative. In deriving it, we equated the GW shrinking time 
to the viscous time of an unperturbed dis c at a radius cor- 
rcsponding to the inner rim of the gap. iTanaka fe Menoul 
(|2010r >. however, showed that the steep density gradient at 
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the inner edge of the disc will shorten the inward diffu- 
sion timescale of the gas, causing a significant delay in the 
disc-binary detachment, therefore increasing the number of 
sources with observable minidiscs. Moreover, in the /3-disc 
picture, the consumption time of a minidisc of ~ 30Rs can 
easily exceed 10 4 years (see equation (jXBJ) , implying that a 
relevant fraction of the detached sys tems may still be signif- 
icantly accreting IjChang et alj|20ich . It is in any case worth 
noticing that the relevant population is sizable (maybe few- 
to-hundred sources) and probably mostly composed by very 
low redshift systems (z < 0.2, see lower left panel of figure 



7.2 Simulations of double Kq line observability 

We assume the iron line to be simil ar to that observed 
from the archetype MCG-6 - 30-15 dTanaka et all 1 19951 ; 
iBrenneman fc Reynolds! l200rj ; iMiniuttil I2007J I We aim to 
assess the feasibility of detecting pairs of relativistic iron 
lines which may be emitted from the inner minidiscs around 
the merging black holes discussed herein, and to use them 
to constrain the properties of the binary system, i.e., black 
holes spin, radial velocity, inclination etc. In order to investi- 
gate the feasibility of utilizing broad Fe Ka emission from a 
binary black hole merg er, spectra were simulated with xspec 
vl2.(jj (JArnaud 1996). We assume the availability of a next 
generation high throughput X-ray telescope, in comparison 
to current instruments, i.e., XMM-Newton. Specifically, we 
use the response matrices created for the proposed Athena 
mission concepil 10 !. The current implementation envisions 2 
focal plane instruments (i) the WFI a wide field CCD imager 



field of view (2.4 arcmin 2 ) calorimeter providing high res- 
olution spectra in the Fe Ka region of 5eV. Hereafter, we 
assume that the binary system has been identified, opening 
the possibility to study it in high spectral resolution with 
the XMS. 

We assume a total mass of ~ 10 9 Mq for the binary 
pair and a mass ratio of 1/3, consistent with the binary 
population presented in Section 17.11 A fiducial redshift of 
0.1 and a bolometric luminosity of 10 % LEdd and an asso- 
ciate d bolometric correction of ~ 20 IjVasudevan fc Fabianl 
l2009l ; lLusscll2010l ). i.e., f x i ~ 10" 11 erg s" 1 cm- 2 ,f x2 ~ 
3 x 10 -11 erg s _1 cm -2 . The parameters of our fiducial sys- 
tem are marked by vertical ticks in figure [9] The black hole 
binary is assumed to have a periodicity consistent with a 
velocity separation of 10 4 km s _1 . For simplicity, we assume 
the energy of the first iron line to be consistent with 6.4 
keV emission at a redshift of 0.1, while the second line is 
blueshifted by 10 4 km s _1 relative to the first line. 

In order to simulate the expected X-ray spectrum from 
such a system, we first make a number of simplifying as- 
sumptions. The broadband continuum emission from each 
accreting black hole is assumed to have a power-law shape, 
where the spectral index has a value of Y — 1.8. The rel- 
ativistic iro n line is modelled using the laor line profile 
(|Laorlll99lr ). where the equivalent width of the line is set 
to be approximately 150eV, as expected for reflection from 



9 http : //heasarc . nasa . gov/xanadu/xspec/ 

10 http : //www . mpe . mpg . de/athena/home . php?lang=en 
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Figure 10. Example of a blind fit model to the simulated spec- 
trum described in i|7,2l The continuum is modeled with a power- 
law (r ~ 1.8). Two relativistic lines in addition to a narrow line 
are required to provide an accurate fit. The second relativistic line 
is required at greater than the 5<r confidence level as measured by 
a simple F-test. The residual panels indicate no lines (bottom), 
a single broad line plus Gaussian (middle) and 2 broad lines plus 
Gaussian. See text for details. 



a disc surrounding a black hole IjGeorge fc Fabianl 1 19911 ) . 
The emissivity profile of the disc is fixed at R -3 , and the 
outer radius of the emitting region is assumed to be 40 
Rs (see £|4.4[) . A narrow line consistent with Fe Ka emis- 
sion from material at much larger radii is also included. 
These values are consistent with current available obser- 
vatio ns (e.g.. lde La Calle Perez et al.ll2010l ; iGuainazzi et al.l 
120111 ) . A more thorough and self-consistent modelling of the 
expected spectrum from a MBHB system, for example, in- 
cluding complex absorbers and self consistently calculating 
the continuum plus r eflected emission is bey ond the scope 
of this work (e.g., see IBrenneman et ajj|2011f ). and as such 
we defer a detailed analysis of this problem for future inves- 
tigations. 

The model described above was defined in xspec 
as pha*((zpo+laor) + (zpo+laor)+ gauss) and 200 spectra 
were simulated in each case, using the latest available re- 
sponse matrices. An example of a simulated spectrum is 
shown in the top panel of figure [101 The interstellar absorp- 
tion component is modelled via the pha model; however, as 
we are interested in the Fe line region, modest values for 
the column density (i.e., < 10 22 cm -2 ) will have a negligi- 
ble effect. Hence the column density was held constant in 
all fits and as such is not discussed further. The spectral 
index F is assumed to be identical in both sources as it is 
inherently difficult to accurately extract both indices cor- 
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Figure 11. Left: Results of the best fit model in figure [TOl for binary system with spins of 0.0 & 0.9 at an inclination of 30° for a random 
sample of 20 model realisations. Right: As in the figure on the left but here the black hole spins are set to 0.3 & 0.7 respectively at 
an inclination of 30°. The colored points indicate the results of fits initialized at the best fit values, whereas the plus & cross symbols 
indicate the results of the best fit blind model. Error bars (90% confidence level) are plotted on the blind model results only for clarity. 
The exposure time is 200ks. In both cases the lines are resolved. 



rectly, if the difference between them is small, due to the 
narrow assumed bandpass (0.5 - 10.0 keV). After each spec- 
trum was simulated the fit & error commands were then 
run in order to measure the best fit parameters given the 
signal-to-noise (SNR) of the spectrum. This results in scat- 
ter around the defined model which is indicated in color in 
figure 1111 This model was then removed and a new 'blind' 
model was defined. In this case only 2 parameters are known 
a priori, the redshift and the line of sight column density. 
Additionally the inclination of the line is fixed at the known 
value. All other parameters are initialized at reasonable val- 
ues given the observed spectrum. The model is first fit with 
only a single relativistic line and the narrow Fe Ka line 
(pha*(zpo+laor+zgauss)). After the best fit is obtained, a 
second relativistic line is added to the model and the signifi- 
cance of this line is calculated using a simple f test, e.g., see 
figure 1101 In all cases the inclination of both lines are tied 
to each other, as we expect the angular momentum of each 
black hole to be aligned by this time in th e merging process 
([Bogdanovic et al J [20071 : iDotti et alj|2010h . 

We focus on 2 primary cases for the spin of the merging 
black holes (i) ai = 0.0, &2 = 0.9, and (ii) ai = 0.3, &2 = 0.7. 
In figure [10] we plot an example best fit model to one of our 
simulated spectra in case (i) above. Here the inclination is 
30°. In figure [TT1 we display the results of the fit to the rel- 
ativistic iron lines. The input model is in colour, whereas 
the subsequent blind fits are indicated by the plus and cross 
symbols. The exposure time in this case is 200ks. It is clear 
that it will be possible to resolve and constrain both Fe 
lines to high accuracy. Similar results are achieved for case 
(ii). At higher inclinations our ability to accurately recover 
the line parameters deteriorates as expected due to the rel- 
ative narrowing of the line profile. We note that in all of 
the lower inclination models (< 40°), all of the blind fits 
require the presence of two relativistic lines at greater than 
the 5a confidence level. As we move to higher inclinations, 
this decreases somewhat whereby at 60° only ~ 95% of our 
realisations require 2 relativistic lines. 



The primary uncertainty in all models is the inclina- 
tion of the binary system, due to the degeneracy in the line 
shape with changing inclination and/or spin. If we relax the 
initial inclination constraint imposed above, our ability to 
accurately determine the line parameters is weakened. For 
example in case (i) above, a second line is required at > 5<r 
level on only ~ 85% of our models, while only 95% require a 
second line at the 3<r level. We expect, however, to get some 
indicative pri or on the system inclinati on from PTA observa- 
tions. Infact. lSesana fe Vecchid (|2010l ) showed that the GW 
signal analysis will enable a determination of the system 
inclination within a ~20 deg accuracy (assuming a source 
SNR of 10). Models were also created with additional narrow 
emission/absorption lines consistent with Fe XXV/XXVI. 
As in the case of the narrow 6.4 keV line, the exquisite reso- 
lution provided by the calorimeter in the Fe K region allows 
these lines to be easily resolved. We also experimented with 
decreasing the velocity separation of the broad line compo- 
nents; however, in this case the results are strongly incli- 
nation dependent. The use of a more accurate relativistic 
line p rofile model, e.g., kerrdisk (JBrenneman fe Reynolds! 
2006), would help in this case. 

There are a number of caveats which should be consid- 
ered when interpreting the results above. For example, even 
in the case of MCG-6-30-15 where long high SNR observa- 
tions exist, there is uncertainty in the paramet ers of the ob- 



served line when modelled by different group s (Fabia n et a.1.1 
120021 : iRevnolds et~aH 120041 : iBrennemaii fe Reynolds! |2006| ). 



In addition, complex rel ativistic effects (e.g., light bending, 
iMiniutti fe Fabianll2004l ) may render the line indistinguish- 
able from the continuum in the absence of highly sophis - 
ticated spectral modelling (e.g.. iBhavani fe Nandra|[201ll ). 
Finally, we note that an alternative physical model has been 
proposed to explain the relativistic lines observed in nu- 
merous AGN. In this scenario, the lines are produced by 
a combination of non-relativistic fluorescent line emission 
from distant material and complex absorption by material 
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in the in ner accretion flow , e.g., for further details see the 
review by I Turner fc Milled IJ2009h . 



proposed strategy for coordinating multimessenger observa- 
tions. We summarize in the following our main results: 



8 DISCUSSION AND CONCLUSIONS 

The distinctive signatures of PTA-MBHBs highlighted 
above open interesting scenarios for combining pulsar tim- 
ing and X-ray observations in the coming years. Even though 
current PTA efforts (PPTA, EPTA, NANOGrav) may even- 
tually succeed in the detection challenge, multimessenger as- 
tronomy prospects are particularly promising in the context 
of the planned SKA. If a nominal Ins sensitivity level will be 
achieved, SKA will make GW astronomy with pulsar timing 
possible. Several hundreds of signals emitted by the low red- 
shift population of MBHBs will add up to form a confusion 
noise, on top of which some (possibly a hundred) sources 
will be individually resolvable and their position in the sky 
determined. 

In this context, X-ray (but also other bands, not consid- 
ered in this paper) observations may play either a prepara- 
tory or a follow-up role. On one hand, the detection of a 
population of periodic X-ray sources in the present decade, 
may provide useful information for PTA observations. By 
the time the SKA will be online, we may have identified 
a catalogue of systems in the sky from which we expect 
to detect GW signals, being also able to give indications 
about the expected frequency. This may substantially facil- 
itate the PTA detection and analysis pipelines, by reducing 
the search parameter space for specific signals. On the other 
hand, several resolvable sources with high SNR may be lo- 
cated in the sky with enough accuracy to make follow-up 
X-ray monitoring possib le. To be conservative, according to 
ISesana fc Vecchid 1)20101 ) . typical PTA source sky location 
accuracy is expected to be in the order of tens of square 
degrees for a source SNR= 10 (but, under specific detec- 
tion condi tions, it can actually be an order of magnitude 
better, see lCorbin fc Cornish! 120101 ). According to figure^ 
typical masses are > few x 10 8 M© at z < 0.3. Such sys- 
tems should be host i n massive galaxies wi th stellar bulges 
of masses > 10 11 Mq IIGiiltekin et alJl2009T ). whose number 
density is < 10~ 3 Mpc~ 3 (|Bell et alj|2003h Considering a sky 
location error of 10 deg 2 , the comoving volume enclosed in 
the error box is 2 x 10 _3 Gpc 3 : maybe up to few thousands 
candidates falls in the error-box. Since the presence of a dis- 
tinctive counterpart relies on the assumption of accretion, 
active galaxies only should be selected, which may leave us 
with few dozens of candidates. Further down selection may 
be performed by keeping galaxies that show signatures of 
recent merger activity [3 If just one or at most a few sys- 
tems are identified, ultra deep X-ray exposure may be used 
to reveal characteristic double relativistic Fe Kq emission 
lines. 

Such investigation will likely be possible in the next 
decade with future generation of pulsar timing arrays (in 
particular with the SKA) and X-ray observatories. In this 
paper we quantified for the first time the population of ex- 
pected sources, and their characteristic signatures, and we 



11 A detailed study of the statist ics of candidate hos ts can be 
found in the independent study bv lTanaka et al.l 112011! ). which is 
complementary to ours in several aspects. 



• About a thousand MBHBs are expected to contribute 
to the PTA detectable signal in the frequency range 10 -9 — 
10 _6 Hz, at a level of Ins or more. Uncertainties in the 
MBHB population model and in the detailed dynamics of 
the contributing systems may impact such figures by a fac- 
tor of a few. 

• We assumed the standard picture of MBHB migration 
in circumbinary discs. For popular geometrically thin, op- 
tically thick discs described by typical parameters, we find 
that many of these PTA sources (order of few hundreds, 
20%-to-50% of the total population, for 0.1 < a < 0.3 and 
0.1 < rh < 0.3) are coupled to their circumbinary discs, 
making a strong case for looking to possible electromagnetic 
signatures rising by the disc-binary dynamical interplay. 

• Detailed SPH simulations of the eccentric binary-disc 
interaction highlight periodic streams of gas leaking from 
the inner edge of the circumbinary disc, through the low 
density central gap. The rate at which such streams feed 
the central binary can be a significant fraction of the Ed- 
dington accretion rate of the MBHB (after initial relaxation 
of the system, we find an average m w 0.4). The stream- 
ing periodicity is extremely sharp, with a variation in the 
MBHB feeding rate of a factor of two-to-four. 

• By modelling analytically the dynamics and the emitted 
spectrum of the inflow- fed minidiscs forming around the two 
MBHs, we identified several observational features that may 
be distinctive of an accreting MBHB. We concentrated on 
the X-ray domain, by separately studying X-ray periodicity 
and double relativistic Fe Ka lines. 

• Assuming all merging MBHBs are accreting, we esti- 
mate about 100-500 (depending on the detailed property of 
the circumbinary disc) periodically variable X-ray sources 
with periods between one and five years. The observed flux 
on Earth for most of the sources is larger than 10 _13 erg 
s _1 cm -2 , within the sensitivity reach of the upcoming X- 
ray all sky monitor eRosita. However, a lightcurve with at 
least 15-20 points (eRosita will point each region of the sky 
at 8 different times only) is required for a statistically sig- 
nificant detection of such a periodicity. 

• Iron line identification requires deep, targeted observa- 
tions. We therefore consider as suitable candidate only those 
PTA sources that are individually resolvable, for which the 
sky location can be identified to some accuracy. The number 
of suitable targets is of the order of few tens, depending both 
on the ability of PTA of resolving sources in the sky, and on 
the details of the disc-MBHB decoupling and the physical 
nature of the minidiscs. 

• Assuming that the accretion flow onto each black hole 
is capable of generating relativistic Fe Ka emission lines, we 
have demonstrated that it will in principle be feasible, via 
high spectral resolution observations with a next generation 
X-ray observatory (e.g. Athena), to identify double Ka line 
features. Such features can be used to estimate the prop- 
erties of the black holes, forecasting, in combination with 
modelling of the PTA signal, the possibility to constrain the 
space-time around the black hole to unprecedented accuracy. 

• Note that we assumed all MBHB to be surrounded by 
a circumbinary disc in their late evolutionary stage, prior to 
merger. Even though this is likely for gas reach merging sys- 
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terns, it might be a too extreme assumption for the massive, 
low redshift sources relevant to our study. Even assuming 
that only 10% of the systems are surrounded by a circumbi- 
nary disc, the number of detectable sources is still interest- 
ing: about 10-50 bright, periodic X-ray sources; at least a 
few individually resolvable PTA sources showing double Fe 
Ka line profiles. 

Even though GW detection of MBHBs remains a challeng- 
ing task for the present and future astrophysical generations, 
systematic pulsar timing campaigns are ongoing, and their 
accuracy and sensitivity will inevitably improve in the com- 
ing years. Pulsar timing will therefore provide a safe, open 
GW window on the low frequency Universe. The combina- 
tion with electromagnetic observations, such the ones pro- 
posed in this paper, will help exploiting GW detection ca- 
pabilities at their best, providing a lot of information about 
the population and dynamics of MBHBs. The present paper 
is just a first step into the realm of multimessenger astron- 
omy with pulsar timing, hoping that investigators from both 
the GW and the X-ray/optical/radio communities will take 
the challenge following our footsteps. 
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